# clear work space
rm(list = ls())

# import packages
library(tidyverse)
library(randtoolbox)
library(maxLik)
library(frontier)

# Halton Draw
R <- 200
halton_sample <- halton(R)

# input rawdata
rawdata <- read_csv("greene2003/TableF4-4.csv")

# model specification
fm <- log(COST/PF)~log(PL/PF)+log(PK/PF)+log(Q)+I(log(Q)^2)

# define estimation function
sfa_gamma <- function(formula, data, ineffDecrease = TRUE, gamma_model = TRUE){
  coefs <- coef(frontier::sfa(formula = formula, data = data, ineffDecrease = ineffDecrease))
  start_value <- c(coefs[1:(length(coefs)-2)],
                   theta = unname(1/sqrt(coefs["gamma"]*coefs["sigmaSq"])),
                   sigmav = unname(sqrt(coefs["sigmaSq"]-coefs["gamma"]*coefs["sigmaSq"])),
                   P = 1)
  
  y <- model.response(model.frame(formula, data))
  X <- model.matrix(formula, data)
  N <- nrow(data)
  
  logL <- function(params){
    betas <- params[1:(length(params)-3)]
    theta <- params["theta"]
    sigmav <- params["sigmav"]
    P <- params["P"]
    
    # residuals
    ep <- y - X %*% betas %>%
      as.vector()
    
    # sample mean
    if (ineffDecrease){
      # production
      mu <- -ep-theta*sigmav^2
    } else {
      # cost
      mu <- ep-theta*sigmav^2
    }
    
    # probability lower bound
    P_L <- pnorm(-mu/sigmav)
    
    # sample draw
    u_sample <- function(mu_i, PL_i){
      mu_i + sigmav*qnorm(PL_i + halton_sample*(1-PL_i))
    }
    
    # sample function average
    msl <- 0
    for (index in seq_along(mu)){
      draws <- u_sample(mu[index], P_L[index])
      msl <- log(mean(draws^{P-1})) + msl
    }
    
    if (ineffDecrease){
      # production
      return(N*P*log(theta) -
               N*log(gamma(P)) + 
               theta*sum(ep) + 
               N/2*theta^2*sigmav^2 + 
               sum(log(1-P_L)) + 
               msl)
    } else {
      # cost
      return(N*P*log(theta) -
               N*log(gamma(P)) -
               theta*sum(ep) + 
               N/2*theta^2*sigmav^2 + 
               sum(log(1-P_L)) + 
               msl)
    }
  }
  
  if (gamma_model){
    return(maxLik(logLik = logL,
                  start = start_value,
                  method = "NM"))
  } else {
    return(maxLik(logLik = logL,
                  start = start_value,
                  method = "NM",
                  fixed = "P"))
  }
}

# do estimation
# gamma model
result <- sfa_gamma(fm, rawdata, ineffDecrease = FALSE)
print(summary(result))

# exponential model
result <- sfa_gamma(fm, rawdata, ineffDecrease = FALSE, gamma_model = FALSE)
print(summary(result))